##########################################################################################################################################################
### Setup

### Source setup
source("code/1_modeling_setup.R")



##########################################################################################################################################################
### Analyze relationship between bicameralism IVs and pr(vote)

### Models

# Run
mods_vote_divided_bicam_int <- vector(mode='list', length=2)
mods_vote_divided_bicam_int[[1]] <- glm(voted_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_vote_divided_bicam_int[[2]] <- glm(voted_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))


### Substantive effect sizes

# DV observed proportion
vote_sd <- bill %>%
  filter(!state %in% high_missingness) %>%
  pull(voted_sles) %>%
  sd()

# Seat ratio
eff_seat_vote <- plot_model(mods_vote_divided_bicam_int[[1]], type='pred', terms=c('seat_ratio', 'divided3'))$data %>% 
  as.data.frame() %>% 
  mutate(x=round(x,1)) %>%
  filter(x %in% c(1.0, 2.4)) %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2.4],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv=abs(eff_size/vote_sd))

# Term difference
eff_term_vote <- plot_model(mods_vote_divided_bicam_int[[2]], type='pred', terms=c('term_diff2', 'divided3'))$data %>% 
  as.data.frame() %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv = abs(eff_size/vote_sd))


### Predicted effects

# Seat ratio
plot_vote_divided_seat <- plot_model(mods_vote_divided_bicam_int[[1]], type='pred', line.size=2, dot.size=5,
                                     terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(x='Seat Ratio', y='Probability of Reaching a Vote', title='', linetype='') + 
  theme_bw() +
  geom_text(data=eff_seat_vote, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.3, vjust=0.5, size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")

# Term difference
plot_vote_divided_term <- plot_model(mods_vote_divided_bicam_int[[2]], type='pred', line.size=2, dot.size=5,
                                     terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Probability of Reaching a Vote', title='', shape='') + 
  geom_line(linetype='dotted', linewidth=1.2, position=position_dodge(width=0.1)) +
  theme_bw() +
  geom_text(data=eff_term_vote, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.2, position=position_dodge(width=0.1), size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")



##########################################################################################################################################################
### Analyze relationship between bicameralism IVs and pr(law)

### Models

# Run 
mods_law_divided_bicam_int <- vector(mode='list', length=2)
mods_law_divided_bicam_int[[1]] <- glm(law_sles ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))
mods_law_divided_bicam_int[[2]] <- glm(law_sles ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(bill, !state %in% high_missingness))


### Substantive effect sizes by party

# DV observed proportion
law_sd <- bill %>%
  filter(!state %in% high_missingness) %>%
  pull(law_sles) %>%
  sd()

# Seat ratio
eff_seat_law <- plot_model(mods_law_divided_bicam_int[[1]], type='pred', terms=c('seat_ratio', 'divided3'))$data %>% 
  as.data.frame() %>% 
  mutate(x=round(x,1)) %>%
  filter(x %in% c(1.0, 2.4)) %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2.4],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv=abs(eff_size/law_sd))

# Term difference
eff_term_law <- plot_model(mods_law_divided_bicam_int[[2]], type='pred', terms=c('term_diff2', 'divided3'))$data %>% 
  as.data.frame() %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==1],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv = abs(eff_size/law_sd))


### Predicted effects 

# Seat ratio
plot_law_divided_seat <- plot_model(mods_law_divided_bicam_int[[1]], type='pred', line.size=2, dot.size=5,
                                    terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(x='Seat Ratio', y='Probability of Becoming Law', title='', linetype='') + 
  theme_bw() +
  geom_text(data=eff_seat_law, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.4, size=6)  +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")

# Term difference
plot_law_divided_term <- plot_model(mods_law_divided_bicam_int[[2]], type='pred', line.size=2, dot.size=5,
                                    terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Probability of Becoming Law', title='', shape='') + 
  geom_line(linetype='dotted', linewidth=1.2, position=position_dodge(width=0.1)) +
  theme_bw() +
  geom_text(data=eff_term_law, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=3.5, size=6)  +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")



##########################################################################################################################################################
### Analyze relationship between bicameralism IVs and proportion of members voting yea

### Models

# Run
mods_yea_divided_bicam_int <- vector(mode='list', length=2)
mods_yea_divided_bicam_int[[1]] <- lm(yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_divided_bicam_int[[2]] <- lm(yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))


### Substantive effect sizes

# DV standard deviation
yea_sd <- votes %>% 
  filter(!state %in% high_missingness) %>% 
  pull(yea_prop) %>% 
  sd()

# Seat ratio
eff_seat_yea <- plot_model(mods_yea_divided_bicam_int[[1]], type='pred', terms=c('seat_ratio', 'divided3'))$data %>% 
  as.data.frame() %>% 
  mutate(x=round(x,1)) %>%
  filter(x %in% c(1.0, 2.4)) %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2.4],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv=abs(eff_size/yea_sd))

# Term difference
eff_term_yea <- plot_model(mods_yea_divided_bicam_int[[2]], type='pred', terms=c('term_diff2', 'divided3'))$data %>% 
  as.data.frame() %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv = abs(eff_size/yea_sd))


### Predicted effects (divided government interaction models)

# Seat ratio
plot_yea_divided_seat <- plot_model(mods_yea_divided_bicam_int[[1]], type='pred', dot.size=5, line.size=2,
                                    terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(x='Seat Ratio', y='Proportion of Members Voting Yea', title='', linetype='') + 
  theme_bw() +
  geom_text(data=eff_seat_yea, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.5, size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")

# Term difference
plot_yea_divided_term <- plot_model(mods_yea_divided_bicam_int[[2]], type='pred', dot.size=5, line.size=2,
                                    terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Proportion of Members Voting Yea', title='', shape='') + 
  geom_line(linetype='dotted', linewidth=1.2, position=position_dodge(width=0.1)) +
  theme_bw() +
  geom_text(data=eff_term_yea, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.2, position=position_dodge(width=0.1), size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")



##########################################################################################################################################################
### Analyze relationship between bicameralism IVs and proportion of minority party members voting yea

### Models

# Run
mods_yea_min_divided_bicam_int <- vector(mode='list', length=2)
mods_yea_min_divided_bicam_int[[1]] <- lm(min_yea_prop ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))
mods_yea_min_divided_bicam_int[[2]] <- lm(min_yea_prop ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), data=filter(votes, !state %in% high_missingness))


### Substantive effect sizes

# DV standard deviation
yea_min_sd <- votes %>% 
  filter(!state %in% high_missingness) %>% 
  pull(min_yea_prop) %>% 
  sd(na.rm=TRUE)

# Seat ratio
eff_seat_yea_min <- plot_model(mods_yea_min_divided_bicam_int[[1]], type='pred', terms=c('seat_ratio', 'divided3'))$data %>% 
  as.data.frame() %>% 
  mutate(x=round(x,1)) %>%
  filter(x %in% c(1.0, 2.4)) %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2.4],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv=abs(eff_size/yea_min_sd))

# Term difference
eff_term_yea_min <- plot_model(mods_yea_min_divided_bicam_int[[2]], type='pred', terms=c('term_diff2', 'divided3'))$data %>% 
  as.data.frame() %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv = abs(eff_size/yea_min_sd))


### Predicted effects

# Seat ratio
plot_yea_min_divided_seat <- plot_model(mods_yea_min_divided_bicam_int[[1]], type='pred', dot.size=5, line.size=2,
                                        terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(x='Seat Ratio', y='Proportion of Minority Members Voting Yea', title='', linetype='') + 
  theme_bw() +
  geom_text(data=eff_seat_yea_min, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), 
                                       group=eff_size_dv), inherit.aes = FALSE, hjust=1.5, size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")

# Term difference
plot_yea_min_divided_term <- plot_model(mods_yea_min_divided_bicam_int[[2]], type='pred', dot.size=5, line.size=2,
                                        terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Proportion of Minority Members Voting Yea', title='', shape='') + 
  geom_line(linetype='dotted', linewidth=1.2, position=position_dodge(width=0.1)) +
  theme_bw() +
  geom_text(data=eff_term_yea_min, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.2, position=position_dodge(width=0.1), size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")



##########################################################################################################################################################
### Analyze relationship between bicameralism IVs and likelihood of party unity votes

### Models

# Run
mods_unity_divided_bicam_int <- vector(mode='list', length=2)
mods_unity_divided_bicam_int[[1]] <- glm(unity_vote ~ seat_ratio*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))
mods_unity_divided_bicam_int[[2]] <- glm(unity_vote ~ term_diff2*divided3 + veto_override_threshold + lineitemveto + maj_sponsor + num_cs + maj_share + pty_diff + mds1 + committee_sponsor + maj_unity + gov_pty_override + as.factor(year), family='binomial', data=filter(votes, !state %in% high_missingness))


### Substantive effect sizes

# DV observed proportion
unity_sd <- votes %>% 
  filter(!state %in% high_missingness) %>% 
  pull(unity_vote) %>% 
  sd(na.rm=TRUE)

# Seat ratio
eff_seat_unity <- plot_model(mods_unity_divided_bicam_int[[1]], type='pred', terms=c('seat_ratio', 'divided3'))$data %>% 
  as.data.frame() %>% 
  mutate(x=round(x,1)) %>%
  filter(x %in% c(1.0, 2.4)) %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2.4],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv=abs(eff_size/unity_sd))

# Term difference
eff_term_unity <- plot_model(mods_unity_divided_bicam_int[[2]], type='pred', terms=c('term_diff2', 'divided3'))$data %>% 
  as.data.frame() %>%
  dplyr::select(x, predicted, group) %>%
  group_by(group) %>%
  summarise(
    eff_size = diff(predicted),
    pred_x = predicted[x==2],
    x_text = max(x)
  ) %>%
  mutate(eff_size_dv = abs(eff_size/unity_sd))


### Predicted effects

# Seat ratio
plot_unity_divided_seat <- plot_model(mods_unity_divided_bicam_int[[1]], type='pred', dot.size=5, line.size=2,
                                      terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(x='Seat Ratio', y='Probability of Party Unity Vote', title='', linetype='') + 
  theme_bw() +
  geom_text(data=eff_seat_unity, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.5, vjust=-0.1, size=6) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")

# Term difference
plot_unity_divided_term <- plot_model(mods_unity_divided_bicam_int[[2]], type='pred', dot.size=5, line.size=2,
                                      terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(x='Difference in Term Length between Chambers', y='Probability of Party Unity Vote', title='', shape='') + 
  geom_line(linetype='dotted', linewidth=1.2, position=position_dodge(width=0.1)) +
  theme_bw() +
  geom_text(data=eff_term_unity, aes(x=x_text, y=pred_x, label=paste0(round(eff_size_dv, 2), ' SDs'), group=eff_size_dv), 
            inherit.aes = FALSE, hjust=1.2, position=position_dodge(width=0.1), size=6)  +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        legend.position="none")



##########################################################################################################################################################
### Combine plots

# Figure 6: Term difference
term_leg <- plot_model(mods_unity_divided_bicam_int[[2]], type='pred', line.size=2,
                       terms=c('term_diff2', 'divided3'), colors='bw') +
  labs(shape='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(shape = guide_legend(override.aes = list(size = 15)))
term_leg_plot <- cowplot::get_legend(term_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))
plot_term <- plot_vote_divided_term + plot_law_divided_term + plot_yea_divided_term + plot_yea_min_divided_term + plot_unity_divided_term + term_leg_plot + plot_layout(ncol=2)
ggsave('results/fig6_term_divided_plot.pdf', width=12, height=15, plot=plot_term)

# Figure 7: Seat Ratio
seat_leg <- plot_model(mods_unity_divided_bicam_int[[1]], type='pred', line.size=2,
                       terms=c('seat_ratio', 'divided3'), colors='bw') +
  labs(linetype='Legend') + theme_bw() +
  theme(legend.text=element_text(size=20), legend.title=element_text(size=20)) +
  guides(linetype = guide_legend(override.aes = list(size = 20)))
seat_leg_plot <- cowplot::get_legend(seat_leg) %>%
  as_ggplot() + theme(legend.key.size=unit(10,'cm'))
plot_seat <- plot_vote_divided_seat + plot_law_divided_seat + plot_yea_divided_seat + plot_yea_min_divided_seat + plot_unity_divided_seat + seat_leg_plot + plot_layout(ncol=2)
ggsave('results/fig7_seat_divided_plot.pdf', width=12, height=15, plot=plot_seat)
